EIGENVALUES FOR TWO-LAG LINEAR DELAY DIFFERENTIAL 

EQUATIONS* 
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Abstract. We consider the class of two-lag linear delay differential equations and develop a 
series expansion to solve for the roots of the nonlinear characteristic equation. Supporting numerical 
results are presented along with application of our method to study the stability of a two-lag model 
from ecology. 
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1. Introduction. Delay Differential Equations (DDEs) naturally arise in many 
fields of science and engineering and accordingly have been the focus of extensive 
study over the years. The classic books by Bellman and Cooke [10J and El'sgol'ts 
and Nor kin [27 describe the significant mathematical results up to the 1960 's, while 
subsequent books by Hale and Verduyn Lunel [33 and Diekmann et al., [26 cover the 
later semigroup-based frameworks. A specific area in which delays arise frequently is 
biological modeling, e.g., due to gestation periods, and there are several good texts 
on population modeling and analysis with DDEs |23l |3"T1 |4"01 |4"3"] including a recent 
introductory textbook by Smith [54 . Lastly, we direct the interested reader to several 
recent review articles [32 , 47] and edited works [5] |32J [53] which thoroughly cover the 
state-of-the-art advances in the study and application of DDEs. 

1.1. DDE Stability Analyses. Consider the first order, linear, single-lag DDE 

(1.1) x(t) = ax(t) +Px(t -r); t > 0, 

(1.2) x(t) = <j>(t);t<Q, 



where (f>(t) is the initial history function (which need not satisfy (l.ll). The Laplace 



transform of (1.1) is 



(1.3) s = a + /3e- ST , 

a transcendental equation with an infinite number of roots {sj} in the complex plane 



(see Figure 3.2 for an example). Thus the solutions are constructed as an infinite 
series 

oo 

y(t) = , 

j=-oo 

and we note that these basis functions e Sjt form a Schauder basis for L p functions 
over a finite domain. The coefficients Cj can be computed via evaluating the history 
function (j){t) and the basis functions e Sjt at 2N + 1 timepoints in [— r, 0] and then 
linearly solving for 2N + 1 of the C/s. 

The nonlinear eigenvalue equation in ( |1.3[ ) belongs to a well-known class of equa- 
tions known as exponential polynomials [5J 1461 148| and many researchers in DDEs 
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1.2 Relationship between DDEs and Lambert W function 1 INTRODUCTION 



have (understandably) chosen to study them [19]. In the broader mathematics 
community, research into this class of algebraic equations is very active and extends 
to numerical schemes for computing the roots (see papers by Jarlebring such as [36 ) 
and applications to fields afar as quantum computing [51J. 

A primary goal in the study of linear (and linearized) DDEs is to identify the 
value of the the delay r for which equilibria become unstable, i.e., when the real part 
of the principal root Re(so) becomes positive. A traditional strategy for proving that 
such a r exists is to separate ( 1.3 1 into real and imaginary equations and prove for the 
t% such that so(t*) is purely imaginary, the eigenvalue passes through to the positive 
real half of the plane, i.e., 



ch 



Re(s (r)) 



> 0. 



Many researchers have been very successful in establishing stability results using this 
framework [30] , even for system with distributed delays [20] [62] , multiple delays |11[ 
STJISU], an d fractional order derivatives |25| . 

1.2. Relationship between DDEs and Lambert W function. In their 2003 
article, Asl and Ulsoy [3] connected the roots of (1.3) with solutions to a specific 
exponential polynomial known as the Lambert W equation |21[ 122], Briefly: to solve 
for the roots to (1.3), rewrite the equation in the form 

t(s - a)e T( - s -^ = p T e- aT , 

and define a new function W such that W(f3Te~ aT ) = t(s — a). The equation to be 
solved is now 



(1.4) 



W(fiTe- aT )e w{fi7 



= (3re- 



which is in the same form as W(x)e w ( x > = x , the well-studied Lambert W equation. 
The roots to (1.3 1 are thus 



1 



W(j3Te- aT ) + a. 



Caratheodory provided the principal solution to the Lambert W function 



W (x) 



OO 

E 

n=l 



and all other branches have been computed as well 
(l.5)W 3 (x) = \n(x) + 2mj - ln(ln(a;) + 2nij) +f^J2 Cl 



ith coefficients 



r(-D' 



1=0 m=l 



/ + m 

l + l 



(ln(x) + 2irikj) 



l+m 



expressed in terms of the unsigned Stirling numbers of the first kind [T]. This ex- 
pansion to compute the different roots for the Lambert W function was presented by 
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de Bruijn [24J, with Comtet [18J rewriting de Bruijn's expansion in terms of the Stir- 
ling numbers and Corless et al. [2.1] correcting a typographic error in equation (2.4.4) 
in Comtet [TH]- Lastly, we note the lambertw function in Matlab, LambertW in 
Maple, and ProductLog in Mathematica allow for the numerically accurate compu- 
tation of Wj(x) for any branch j. 

Since Asl and Ulsoy's original article, the connection between DDE eigenvalues 
and solutions to Lambert W has been substantially extended and expanded |13( HH 
HS1 LTSl HTJ EH [351 EZl [52J [551 [53 EH [601 El], culminating in a monograph [59] . These 
results, however, do not extend to multi-lag DDEs as (to the best of our knowledge) 
the nonlinear eigenvalue problem can be cast as a Lambert W equation only when 
the DDE h single lag. 

In this work, we develop a series expansion for the roots to the nonlinear eigenvalue 
problem corresponding to a two-la,g DDE. While we were inspired by the original 



derivation of (1.5), the development presented below is a bit more involved than 
the Lambert W derivation in de Bruijn [24 partially due to our goal of making the 
truncated sums computationally efficient to evaluate. 

Section [2] presents the main result of the paper, deriving the formula for the 
eigenvalues. Section [3] provides numerical results supporting our claims of accuracy 
and convergence and including an application to a two-lag model from ecology. Lastly, 
Section [4] summarizes the results of the work and discusses the natural next steps for 
future work. 

2. Two-Lag Development. In this section, we develop our series expansion for 
the roots to the nonlinear characteristic equation arising from the two-lag DDE. Note 
that for the multi-delay case, the conventional approach focuses on developing bounds 
for the roots (see [10 , §12] for a good summary and overview), whereas here, we are 
explicitly computing the roots. 

Consider a DDE with two lags n and T2 

x(t) = ax(t) + f3x{t — n) + jx(t — r 2 ) , 

with corresponding coefficients a, /3, and 7. We will rescale and rearrange terms, 
defining several new variables and parameters in order to reduce the equation to 
its essential mathematical features. We rescale time by the first lag t = t/ri, the 
parameters so that ot\ = ar%, fix = /3tx, 71 = 7x1 , r = t i/t 2 , and denote x'(t) = 
4fX(i) = Tx4:x{t) to obtain the equation 

x'(t) = axx(t) + (3ix(i - 1) + 7i.x(t - r) . 

After a Laplace transform, the corresponding nonlinear characteristic equation for 
s G C is thus 

(2.1) s = at +Pie- S +j lS - ST . 

We then let A = s — aj., P2 = /?ie _ai , 72 = 7ie~ Qir , and multiply through by e Ar , 
yielding the much cleaner equation 

(2.2) Ae AT =/3 2e ^- 1 ) A + 72 . 

Inspired by the original expansion in [24J, we note that for I72I near to or far from 
zero, we can write A as a small perturbation u of In 72 

A = (In 72 +u)/t, 
3 
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2 TWO-LAG DEVELOPMENT 



and then safely assume that 
(2.3) 



M < II1172I 



Substituting this form of A into equation (|2.2| we find that 
1 



(2.4) 



1 



72 e 



r \ In 72, 
To proceed, we will also assume that 



— } (e ")(^) + 
In 72 In 72 



(2.5) 



1/t, 
72 m 72 



< 1, 



in addition to (2.3). This restriction on the coefficients is a result of the chosen 



rearrangement of (2.1) into (2.2). There is nothing particularly special about (2.5 1 



and choosing to recast (2.1 1 differently would simply necessitate different assumptions 



on the parameters in order for an expansion to converge. In our experience, this is a 
reasonable assumption encountered in several models (see the examples in Section [3]) . 
Moving forward, we can solve (2.4 1 for u as 



u ~ In r + ln(Yin7 2 ) 

which means that 

(2.6) u = lnr + ln(Vin 72 ) +v 

where our attention now focuses on solving for a new variable v. As discussed by 



Corless et al. [ST], the nested logarithms in (2.6 1 need not select the same branch. As 
Corless et al. do, we will let the outer logarithm be the principal branch and denote 
Ln as the inner logarithm for which we have not yet chosen a branch, i.e., 

u = lnr + In( 1 /Lii7 2 ) + v . 



By substituting u back into ( 2.4 ) , we can cancel several terms and simplify our problem 

(^) 



to one of finding the roots v to the equation 

1 



1 



lnr ^ In(!/Ln7 2 ) 



L1172 



Ln72 



L1172 



P2 
72 



72 T 

Ln-,2 



If we let 



1 



Ln.72 ' 



{ 72 r 
\Ln72 



(^) 



lnr In( 1 /Ln7 2 ) 
Ln7 2 



c, 



then our efforts will focus on identifying the roots of 
(2.7) f(v) ee e" v + ce~ v/T - av - 1 - 



Note that the motivation for adding and subtracting c will become self-evident in the 
proof for the following lemma. 

Lemma 2.1. If there exists positive numbers a and b such that \a\ < a and < a 
then (2.1) has a single root inside the region \v\ < b. 
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Proof. Let S denote the lower bound of |e v — 1| on \v\ = b and 5 T denote the 
lower bound of \ce~ v / T — c| on \v\ = b. We can arbitrarily choose b = 7rmin{l,r} 
since e~ v — 1 has only one root at zero for b g (0, 2tt) and c(e~ v ^ T — 1) also has 
only one root at v = for b g (0,27rr). Thus there is only one root at v — for 
e~ v — 1 + c[eT v l T — 1) in the region \v\ < 7rmin{l,r} and that on |u| = 7rmin{l,T} 
we have a lower bound of min{i5, S T }. 

Let a = min{<5, S t }/(2(tt + 1)) and then for |<r| < a, < a, and = 7rmin{l,T} 
we have that 

, i ^ „ min{<5, S T ] min{<5, <5 r } min-M, (5 T | . rr 

\av + m| < |a| H + H < -^-^ • * + = < mm{6, S T } . 

Since \e~ v — \ + ce~ v l T — c\ > min{5, <5 T } on the curve defined by |u| = tt min{l, r}, 
we have satisfied the hypotheses of Rouche's Theorem [12 and can therefore conclude 
that there is only one root to (2.7 1 inside \v\ = 7rmin{l,r}. □ 

By Cauchy's theorem, we can then compute v using the well-known formula from 
complex analysis 

( 8) " 2m f lcl=b /(C) dC " 2m f Kl=b e -C + ce -C/r — a( — 1 — c — fx ' 

for b = 7rmin{l,r}. The 1//(C) can be expanded into the following (absolutely and 
uniformly convergent) power series 

1 OO CO 

( 2J ) 7777 = E E ( e ~ C + ce ~ C/T - 1 - c)- fc - m - 1 C^V m a™ +fc , 

fe=0m=0 



with some coefficients a m . The substitution of (2.9 1 into (2.8 1 and evaluation of the 
contour integral generates a double series in m and k. Note, however, that the m = 
term is zero because the integrand has the same number of roots (at £ = 0) in the 
numerator and the denominator. The doubly infinite and absolutely convergent power 
series for v in a and /x is thus 

oo oo 

(2.10) W = ^^a fcm( rV\ 

fe=0m=l 

where the ak m are coefficients independent of a and The coefficients a^m are com- 
puted by first computing a power series expansion in a and then applying the Lagrange 
Inversion Theorem [1_ to obtain v as a function of fi. To facilitate rearrangement of 
this series into the form of (2.10), we denote 

/ 2 H = /H-/(0), 



and note that the derivative of f(w) with respect to a is — w. Equation (2.101 can 
thus be written (after the power series expansion in a around zero) as 



(2.H) v = j2J2 m * hm > 



\i a 

' ml kl 

k=0 m=l 
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2 TWO-LAG DEVELOPMENT 



where m k is a rising factoriaj^] and h m ,k — lhn„,_>o h mt k(w) with 

Am— 1 . . 

Straightforward application of the product rule yields that 



, l J Kdw" 1 - 1 - 1 J \dw l 
1=0 x 7 v 7 v 

The /th derivative of /2(u>)~' m+fc ) is computed using Faa di Bruno's formula (38] for 
higher derivatives of composite functions and thus 

m—X I 
1=0 p=0 

where 

a k im = i j J ^ + i + i)! ' bfcmp = (m + k) ' 

and B/ p are the partial Bell polynomials (21 IH1 El Recall that for an arbitrary 
sequence {x{\ and (,pgN, the partial Bell polynomials are 

B i , p (x 1 ,..,, ; _ p+1 ) = E Jl!j2! ... J; _ p+l! ( T! j UJ •••((T^TDiJ ' 

with the sum taken over all sequences {ji, . . . ,j n -k+i} G N such that h = m 

and YH=i +1 = k (& simple Diophantine linear system). The sets of indices {j^} 
satisfying these sums are just a way of describing all possible partitions of a set of 
size k [28_. 

In what follows, the argument for the Bell polynomial frequently takes the form 
{/a To improve the clarity of the formulas, we denote 

F„ = {/ 2 (1+i) (0)}? =1 ;nGN, 

to be a finite sequence of higher derivatives of fa. 

Formally taking the limit for lim^^o h m> k{w) yields 



(2.1* 



= O'klva 2~2 p =o l\liq=0 ^klmpq^2m—2—l-q,m—l-p 



where 

(2m + k -IV, , ,,. 
aklm = a>klm\ n , „ K + l + ll 
V 2m — 1 — 21 



^klmpq — bkmp 



/2to - 2 - A (to- 1 -p) 



9 /(2m + fc-l)!' 



2 A rising factorial m fc is defined as m(m + 1) • • • (m + k — 1) and a falling factorial m— is defined 
as m(m — 1) ■ ■ ■ (m — k + 1). We note that in the combinatorics literature, m' 1 ) is more commonly 
used as a rising factorial (also known as the Pochhammer symbol) and (m)j, as a falling factorial. 
Here we choose the notation "m fc " for rising factorial to prevent confusion later on when using the 
superscript/parentheses notation for differentiation. If it also the notation advocated by Knuth |39l . 
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The evaluation of most terms in (2.121 is computationally str aight forward and 
efficient, except the term involving b|^. For example, the sum in (2.11 1 for rn = 10 
and k = 100 can take around one minute on an Intel i7 cpu when using Mathematica 
to analytically find the qth derivative of B; p and then sum the terms. In |45| , Mishkov 

derived a formula for b[^, but to implement it would consists of solving a cascade of 
linear Diophantine equation^ to obtain the correct indices for the sums. 

We will use a recursive property of fi| ^ to avoid having to solve systems of 
Diophantine equations. Note that for an arbitrary sequence X = {xj} 



d 



(X) 



n-j, P -i 



(X), 



and thus in the context of (2.121 we have that 



(2.13) 



n=l V ' r 2 =0 V - / 



recursively allows (eventual) computation of the derivatives by evaluation of the par- 
tial Bell polynomials themselves]^] The use of Mathematica's memoization features 
for recursive sums substantially sped up the computation as well, naturally at the 
expense of increased memory usage 



(2.14) 
where 



We now have that the principal branch of the roots of (2.2 1 is 

1 



A = (In 72 +lnr + ln(- ) + v)/t. 

In 72 



00 00 ,,m „k 

(2.15) v = E E mkh 



' ml k\ 

fc=0m=l 



We note that for i > 1 

/ 2 (1+l) (o) = (-ir(i + c o. 

and thus h rn k is defined as 



i = Q Lp=0 L F akhnpq»2m-2-l-q,rn-l-p(r m -l-q+p)»i y p{J"l-p+l) 
(-(l + CT- 1 )) 2 ™^- 1 

where 

F n = {(-ir(l + ar- n m =1 , 



3 These have no closed form solution and would thus also need to be solved. 

4 For a good summary of known properties of Bell polynomials, we direct the interested reader 
to Aldrovandi [2] §13 & Appendix A. 5]. 

5 Memoization is a technique for reducing the number of function calls by storing the results of 
previous calls 1441 . 
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3 NUMERICAL RESULTS 



(— l) p (k + m)r(m)r(m — p)T(k + m + p) 
qll\T(2 + k + l)T(-l + m)r(-l - I + 2m - q) 



L1172 ' 



( 72 T 
\Ln72 



(^) 



In r + In( 1 /Ln7 2 ) 
L1172 



and recalling that is defined recursively in (|2.13l 



The full description of the jth branch root for the characteristic equation in (2.1 1 
is therefore 



-(ln i 72+ln(— — . 
t Irij 72 



lnr 



0L\ 



where lnj(z) = \z\ + arg(z) + 27rij, v is defined in (2.151, and the subscript on the v 



indicates the chosen branch for the logarithms in a, c, and fj,. 

For the convenience of the reader, we have included the Mathematica code in 
Appendix [A| which solves for sj as a function of a, f3, 7, ti, T2, j, m, and fc. Lastly, 
upon request, the author will supply versions of the code with arguments a.\, /?2, 
72, t, j, m, and k as well as versions which obtain significant speedups by using 
Mathematica's memoization, compilation, and parallelization features. 

3. Numerical Results. In what follows, we denote §j as the actual root. Where 
relevant, the Sj values were computed using the damped Newton's method imple- 



mented in Mathematica's FindRoot applied to equation (2.1 1 



3.1. Convergence. In this section, we provide numerical results for the estima- 
tion of the eigenvalues. Consider the two lag DDE 



(3.1) 



x(t) = ax(t) + Px{t - Ti) + jx(t -T 2 );t>0, 
x(t) = 4>{t) ; t < 0, 



system from above with transfer function 

(3.2) [s-a- (5e- STl - 1 e- ST >)- 1 . 

We note that the accuracy of our asymptotic expansion rests on the assumption that 



(2.5) is small and so we choose 

/3 2 = 0.001,72 = 6,t = 4, ai = -1 



(3.3) 



which means that the parameters in (3.1l are 

(3.4) a = -1 , ; P w 3.68 x 1CT 4 , 7 = 0.1989 , n = 1 , r 2 = 5 , 

and thus 

- w 3.6 x 1CT 4 . 



Vt 1 

72 m 72 



Numerical investigations suggested that the series converges in m about a factor of four 
times faster than it converges in k. Accordingly, Figure |3.1| depicts the convergence 



in m with k 



We note that up to m 



the series converges nicely. For m > 8, 
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3.2 Example 



-2 



-4 



10 



Figure 3.1. Plot o/log 1() |so — sol for the series approxima tion to so as a function of m where 
k = m 4 . The parameters for the two-lag DDE were those from \3.4\ ■ 



however, the series exhibits a reduction of accuracy, likely due to the numerical issues 
inherent in evaluating such large series. 

Lastly, in subsequent computations of Sj, we choose to truncate the series at 
to = 8 and k — 1000 to maintain accuracy. We obtain the following roots as depicted 
in Figure |3.2| For comparison, Figure |3.3| depicts a contour plot of the transfer 
function in (3.2 1 (lighter colors are larger values) and we note that our computed 



roots match up nicely with the singularities of the transfer function. 

3.2. Example. In modeling the life cycle of the Australian blowfly Lucila cup- 
rina, Braddock and van den Driessche [TT] developed the following two-lag DDE 



(3.5) 



dx 
~dt 



= rx(t)(l - a\x{t - n) - a 2 x(t - t 2 )) 



The positive equilibrium is at x* = l/(ai + 02) and a linearization about this point 
(shifting to y = x — x*) yields the equation 



(3.6) 



dt 



-rx*aiy(t - n) + -rx*a 2 y(t - r 2 ) . 



Theorem 6 in Ruan (33] summarizes the results of [TT] nicely, providing 3 different 
sufficiency conditions for the existence of a value of t 2 which will incur a Hopf bi- 
furcation. We present here only the last of the possible sufficiency conditions in the 
following theorem from [IS], since these restrictions also satisfy our conditions for 
accurate series expansion. 

Theorem 3.1. // a\ = a 2 a nd tj > (2x*ra\ 
when t 2 = T2 the two-lag DDE (3.5) exhibits a Hopf bifurcation at the equilibrium 
x* = l/(ai + a 2 ). 

To illustrate the match between our expansion and this theory, let r = 1, a\ = 
a,2 = 1, T\ = 10. Figure [3T4] depicts the principal root as a function of r 2 and we note 
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then there is a > such that 
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-0.9 



-0.8 



-0.7 



-0.6 



-0.5 



-0.4 



-0.3 



-4 - 



Figure 3.2. Eigenvalues Sj computed using the expansion in Section The parameters for 
the two-lag DDE were those from (3.4]) ■ 




Figure 3.3. Contours of log 10 \(s — a — j3e 



■ 7e 



~ ST2 ) 1 | f or c £ C. The horizontal axis 



is the real part of s and the vertical axis is the imaginary part of s. Note how the peaks (lighter 
regions) coincide with the computed roots in Figure \3. S\ The parameters for the two-lag DDE were 
those from 



the root at T2 = 0.379414. We also note that the relation which must be small (2.5 1 
is on the order of 10 -17 for this set of parameters. Figure 3.5 depicts asymptotically 
stable oscillations (solid curve for — 0.2) and unstable oscillations (dashed curve 
for T2 = 1.3). To better illustrate the stable and unstable behaviors, we arbitrarily 
choose an initial history of <fi(t) = 1 ; t < for the stable solution and an initial history 
of <j)(t) — !/io ; t < for the unstable solution. 
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4. Conclusion and Discussion. In this work we have developed an asymptotic 
expansion for the roots to a nonlinear eigenvalue problem associated with two-lag 
DDEs. We have provided numerical evidence supporting the accuracy of our expan- 
sion. We have also applied the expansion to a two-lag DDE model of the Australian 
blowfly, illustrating the Hopf bifurcation transition from stable to unstable oscilla- 
tions. 

There are many intriguing directions one could take with this work. Concerning 
systems of single-lag DDEs, it is known that the approach presented in Asl and 
Ulsoy's original article [3] is only valid when the coefficient matrices in front of the 
state and lagged state commute. This limitation has been discussed in the literature 
[4jl63] with Jarlebring [37. also noting that the matrices need only be triangularizable. 
Fortunately, there is an algorithm for working around this problem and the monograph 
[55] has a very clear description of how to solve single-lag DDE systems using a matrix 
Lambert W function. And a reasonable extension would be to develop the framework 
to identify eigenvalues for two-lag systems of DDEs. Indeed, we have made preliminary 
efforts in this direction as it is a natural extension to some of the author's previous 
work [7]. 

Another direction for future work would be to consider more general multi-lag 
DDEs. We were surprised to discover that the restrictions placed on the coefficients 
in the two-lag DDE were satisfied for many of the DDEs we considered. We therefore 
plan to investigate if the coefficient restrictions for an asymptotic expansion in the 
three-lag (and higher) DDEs are as generous. 

5. Acknowledgements. The author wishes to thank Andrew Christlieb, Anna 
Gilbert, and Peter Petre for the original discussion which led to the motivating ques- 
tion for this paper. Portions of this work were performed with support from DARPA, 
AFRL, and HRL, under contract FA8650-11-C-7158. 
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6 DISCLAIMER 




Figure 3.5. Plot of stable (solid curve for T2 
solutions to the linearized two-lag blowfly model in 
is (f>(t) = 1 for t < and for the unstable solution, 



= 0.2) and unstable (dashed curve for T2 = 1.3) 
J[3. For the stable solution, the initial history 
the initial history is 4>(t) = 0.1 for t < 0. 



6. Disclaimer. The views expressed are those of the author and do not reflect 
the official policy or position of the Department of Defense or the U.S. Government. 
This is in accordance with DoDI 5230.29, January 8, 2009. 

Appendix A. Mathematica Code. 

Here we present the Mathematica code which was used to generate the roots Sj : 
(*DM Bortz*) 

(*dmbortz@colorado . edu*) 
(^Mathematical Biology Group*) 
(*http: //mathbio . Colorado . edu*) 
(*Applied Mathematics Department*) 
(*University of Colorado*) 
(*Boulder, CO 80309-0526*) 
(*June 2012*) 

Clear [Logn, f2p, z, j, w, taul , tau2, alpha, beta, gamma, k, 1, m, n] 
Clear[p, q, mmax, kmax, DBellYq, sigma, c, mu, sj , sjMem] 
Logn[z_,j_] := Log[Abs[z]] + Arg[z] + 2*Pi*I*j ; 
f2p[w_, taul_, tau2_, j_, c_] := (-l)~j*(E~(-w) + 

c* (1/ (E~ (w* (taul/tau2) ) * (tau2/taul) ~j ) ) ) ; 
DBellYq [1_, p_ , q_ , taul_, tau2_ , c_] := 

If [1 < I I p < I I (p == && 1 > 0) I I 1 < p, 0, 
If [q == 0, 

BellY[l, p, Array [f2p[0, taul, tau2, #1, c] &, {1-p+l}] ] , 
Sum [Binomial [1 , r] * 

Sum [Binomial [q - 1, s] * 

12 

Approved for Public Release. Distribution Unlimited. 



REFERENCES 



REFERENCES 



DBellYq[l-r,p-l , q-l-s , taiil ,tau2,c] * 
f2p[0, taul, tau2, 1 + r + s, c] , 
{s, 0, q - 1}] , 
{r, 1, 1 - p + 1}] 

] 

]; 

sigma [alpha_ , gamma_ , taul_ , tau2_ , j _] : = 

l/Logn[gamma*taul*Exp [(-alpha) *tau2] , j] ; 
c [alpha_ , beta_, gamma_ , taul_, tau2_, j_] := 
(beta* (Exp [alpha* (tau2 - taul)] /gamma) ) * 
(gamma*tau2* (Exp [( -alpha) *tau2] / 

Logn[gamma*taul*Exp[(-alpha)*tau2] , j]))~(l - taul/tau2) ; 
mu[alpha_, beta_, gamma_ , taul_, tau2_ , j_] := 

(Log[tau2/taul] + Log [1/Logn [gamma*taul*Exp [ (-alpha) *tau2] , j]]) / 
Logn [gamma*taul*Exp[( -alpha) *tau2] , j] - 
c [alpha, beta, gamma, taul, tau2, j] ; 
s j [alpha_ , beta_, gamma_ , taul_, tau2_ , j_, mmax_ , kmax_] := 
(Logn [gamma*taul*Exp[ (-alpha) *tau2] , j] + Log[tau2/taul] - 
Log [Logn [gamma*taul*Exp[ (-alpha) *tau2] , j]] + 
Sum [Sum [ (Sum [Sum [Sum [ 

(-l)~p Pochhammer [m,k] *( (k + m)*Gamma[m] *Gamma[m-p] * 
Gamma [k+m+p] ) / (1 ! *q! *Gamma [2+k+l] *Gamma [-1+m] * 
Gamma [-1 - 1 + 2*m - q] ) * 
BellY[2*m -2-1-q, m-l-p, 

Array [f 2p [0 , taul , tau2 , #1 , c [alpha, beta, gamma, taul , tau2, j] ]& , 
{Max[0, m-l-q + p]}]]* 

DBellYq[l ,p , q, taul ,tau2 , c [alpha, beta, gamma, taul ,tau2, j] ] , 
{q, 0, 2*m - 2 - 1}], {p, 0, 1}], {1, 0, m - 1}] / 
(f2p[0, taul, tau2, 1, 

c[alpha, beta, gamma, taul, tau2, j]]~(2*m + k - 1))) * 
mu [alpha, beta, gamma, taul, tau2, j]~m * 
(sigma[alpha, gamma, taul, tau2, j]~k/m!/k!), 
{m, 1, mmax}] , {k, 0, kmax}] ) / (tau2/taul) + 
alpha*taul ; 
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